setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
library(openxlsx)
file1=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",colNames = TRUE,sheet=6,startRow = 2)
file2=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",colNames = TRUE,sheet=7,startRow = 2)
file3=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",colNames = TRUE,sheet=8,startRow = 2)
file4=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",colNames = TRUE,sheet=9,startRow = 2)
file5=read.xlsx("E:/0 公共数据库差异情况/db_for_5hmc/1.xlsx",colNames = TRUE,sheet=10,startRow = 2)
file1$unitID=paste(file1$Chr.x,file1$Start.x,sep=":")
file1$source.group="ASoC.iPS"
file2$unitID=paste(file2$Chr.x,file2$Start.x,sep=":")
file2$source.group="ASoC.NPC"
file3$unitID=paste(file3$Chr.x,file3$Start.x,sep=":")
file3$source.group="ASoC.iN_Glu"
file4$unitID=paste(file4$Chr.x,file4$Start.x,sep=":")
file4$source.group="ASoC.iN_GA"
file5$unitID=paste(file5$Chr.x,file5$Start.x,sep=":")
file5$source.group="ASoC.iN_DN"
asoc=rbind(file1[,20:21],file2[,20:21],file3[,20:21],file4[,20:21],file5[,20:21])

AJHG=read.table("E:/0 公共数据库差异情况/db_for_5hmc/2704_ASMs_hg18ToHg19_AJHG.bed",head=F,sep="\t")
AJHG$unitID=paste(AJHG$V1,AJHG$V2,sep=":")
AJHG$source.group="AJHG"

### 22 loci
id1=read.table("../20201104_at.least.one.is.AShM.beta0.pic/normal_down_tumor_down.intersect.txt",head=T,sep="\t")
id2=read.table("../20201104_at.least.one.is.AShM.beta0.pic/normal_up_tumor_up.intersect.txt",head=T,sep="\t")
id1=tidyr::separate(id1,unitID,into=c("chr","position","ref","alt"),sep=":")
id2=tidyr::separate(id2,unitID,into=c("chr","position","ref","alt"),sep=":")
id1=paste(id1$chr,id1$position,sep = ":")
id2=paste(id2$chr,id2$position,sep = ":")
loci22=unique(c(id1,id2))
### 252 loci
rt=read.table("../20201110找到有差异的且具有相同pattern的位点/DC.BF1.intersect.txt",head=T,sep="\t")
rt$num=rowSums(rt[,seq(5,25,4)]>=0,na.rm = TRUE)
rt=rt[rt$num>1,]
rt=tidyr::separate(rt,unitID,into=c("chr","position","ref","alt"),sep=":")
loci252=paste(rt$chr,rt$position,sep=":")

test=asoc[as.character(asoc$unitID) %in% loci22,]#没交集
test=asoc[as.character(asoc$unitID) %in% loci252,]
test=AJHG[as.character(AJHG$unitID) %in% loci252,]
test=AJHG[as.character(AJHG$unitID) %in% loci22,]
